home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / VideoToolboxSources / Binomial.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-07  |  8.0 KB  |  275 lines  |  [TEXT/KAHL]
  1. /*
  2. Binomial.c
  3. Copyright (c) 1990,1991,1992 Denis G. Pelli
  4.  
  5. Various routines that deal with Binomial statistics, including generating random
  6. samples and computing confidence intervals.
  7.  
  8. This file contains two generations of routines that generate confidence
  9. intervals. The original routines BinomialUpperBound and BinomialLowerBound are
  10. old, but are more robust. The new routines are based on a Numerical Recipes in C
  11. routine to compute the incomplete beta function. The new routines add little to
  12. the old routines, and require the Numerical Recipes in C, so I commented them
  13. out. The new routines do have the advantage of allowing an arbitrary confidence,
  14. whereas the old routines were hard coded to produce a 95% confidence interval.
  15.  
  16. HISTORY:
  17. 1/5/91 Added new binomial routines. These routines may be 
  18.         commented out by setting NEW_BINOMIAL to zero.
  19. 1/7/91 Rewrote BinomialUpperBound & BinomialLowerBound to allow specification of the
  20.         desired confidence of the interval.
  21. 4/24/92    dgp    Added BinomialSampleQuickly().
  22. 4/27/92    dgp    Oops. I just noticed that this file was always including nr.h, which you'd
  23.             only have if you own the Numerical Recipes. I moved the include statement
  24.             down inside the conditional below, where it belongs. 
  25. */
  26. #include "VideoToolbox.h"
  27. #include <math.h>
  28.  
  29. #define NEW_BINOMIAL 0
  30.  
  31. long BinomialSample(double p,long n)
  32. /*
  33. Returns a a sample from a binomial distribution: number of heads in n flips of
  34. coin with a probability p of heads on each flip. This is fine if n is small. However,
  35. if n is large, e.g. n>100, then this routine will be rather slow, and you may prefer
  36. to use the more elaborate Numerical Recipes bnldev() routine.
  37. */
  38. {
  39.     long k,i;
  40.  
  41.     k=0;
  42.     for(i=0;i<n;i++) if(p>UniformSample())k++;
  43.     return k;
  44. }
  45.  
  46. int BinomialSampleQuickly(int n)
  47. // Does n coin flips and returns the number of heads. Very fast.
  48. {
  49.     register int i,k=0;
  50.     register short r;
  51.     
  52.     for(i=n;i>=8;i-=8){
  53.         r=randU();    // Use only the upper 8 bits, which are reputed to be more random.
  54.         if(r<0)k++;
  55.         r<<=1;
  56.         if(r<0)k++;
  57.         r<<=1;
  58.         if(r<0)k++;
  59.         r<<=1;
  60.         if(r<0)k++;
  61.         r<<=1;
  62.         if(r<0)k++;
  63.         r<<=1;
  64.         if(r<0)k++;
  65.         r<<=1;
  66.         if(r<0)k++;
  67.         r<<=1;
  68.         if(r<0)k++;
  69.         r<<=1;
  70.     }
  71.     if(i>0){
  72.         r=randU();
  73.         for(;i>0;i--){
  74.             if(r<0)k++;
  75.             r<<=1;
  76.         }
  77.     }
  78.     return k;
  79. }
  80.  
  81. /*
  82. BinomialUpperBound and BinomialLowerBound return a P confidence interval for the
  83. underlying binomial probability assumed to have generated the data.
  84.  
  85. The formula is based on a Gaussian approximation, solving for the p that will
  86. put the observed result the right number of standard deviations away from the
  87. mean, plus or minus 0.5, as a "continuity correction".
  88.  
  89. This is the best formula I could find in looking through several statistics
  90. books. However, the result has to be taken with a grain of salt because it is
  91. not possible to produce a binomial confidence interval that will satisfy the
  92. strict definition of a confidence interval, namely one that will have the
  93. specified probability P of containing the unknown but fixed parameter p. That's
  94. because, unlike the Normal distribution, the Binomial distribution is not
  95. translation invariant. E.g. if p is actually zero then the true confidence interval
  96. ALWAYS contains p.
  97.  
  98. I gave some thought to taking a Bayesian approach, assuming a uniform prior pdf
  99. for p and then computing an a posteriori confidence interval. This can be done,
  100. but it's hard. I concluded that it's pointless because the uniform pdf
  101. assumption is usually unwarranted.
  102.  
  103. Even though this confidence interval is unsatisfactory (i.e. fundamentally
  104. false) it is useful in practice since it typically behaves similarly to the
  105. Normal case, which is theoretically sound.
  106. */
  107. double BinomialLowerBound(double P,long k,long n)
  108. /*
  109. Arguments are the confidence P, and the number k of heads in n flips of a coin.
  110. The returned value p is the lower end of a P confidence interval on the
  111. underlying probability of a head on a single trial.
  112. */
  113. {
  114.     double right,s,ss,p;
  115.     
  116.     if(k>0 && k<n){
  117.         right=k-0.5;
  118.         s=InverseNormal(sqrt(P));
  119.         ss=s*s;
  120.         p=(right+0.5*ss-s*sqrt(right*(1.0-right/n)+0.25*ss))/(n+ss);
  121.         return p;
  122.     }
  123.     if(k==0) return 0.0;
  124.     if(k==n) return pow(1.0-P,1.0/n);
  125.     return sqrt(-1.0);    /* domain error */
  126. }
  127.  
  128. double BinomialUpperBound(double P,long k,long n)
  129. /*
  130. Arguments are the confidence P, and the number k of heads in n flips of a coin.
  131. The returned value p is the upper end of a P confidence interval on the
  132. underlying probability of a head on a single trial.
  133. */
  134. {
  135.     double right,s,ss,p;
  136.     
  137.     if(k>0 && k<n){
  138.         right=k+0.5;
  139.         s=InverseNormal(sqrt(P));
  140.         ss=s*s;
  141.         p=(right+0.5*ss+s*sqrt(right*(1.0-right/n)+0.25*ss))/(n+ss);
  142.         return p;
  143.     }
  144.     if(k==0) return 1.0-pow(1.0-P,1.0/n);
  145.     if(k==n) return 1.0;
  146.     return sqrt(-1.0);    /* domain error */
  147. }
  148.  
  149. #if NEW_BINOMIAL    /* new routines that require Numerical Recipes in C */
  150.     #include <nr.h>                /* prototype for betai() */
  151.     #if 0
  152.         void main(void)
  153.         /* a quick and dirty driver to test some of these routines */
  154.         {
  155.             double p,pUpper,pLower,P,PUpper,PLower;
  156.             int i,n,k;
  157.             
  158.             Require(0);
  159.             n=10;
  160.             p=.2;
  161.             pLower=InverseBinomial(0.05,21,50);
  162.             printf("%f, Abramowitz & Stegun p. 960 Example 18 say this should be 0.3003.\n",pLower);
  163.             P=0.95;
  164.             for(i=0;i<10;i++){
  165.                 k=BinomialSample( p, n);
  166.                 pLower=BinomialLowerBound(P,k,n);
  167.                 pUpper=BinomialUpperBound(P,k,n);
  168.                 PLower=1.0-Binomial(pLower,k,n);
  169.                 PUpper=Binomial(pUpper,k+1,n);
  170.                 printf("n %4d p %6.4f %6.4f≤%6.4f≤%6.4f  %6.4f< <%6.4f\n"
  171.                     ,n,p,pLower,k/(double)n,pUpper,PLower,PUpper);
  172.                 pLower=InverseBinomial(1.0-sqrt(P),k,n);
  173.                 pUpper=InverseBinomial(sqrt(P),k+1,n);
  174.                 PLower=1.0-Binomial(pLower,k,n);
  175.                 PUpper=Binomial(pUpper,k+1,n);
  176.                 printf("  %4d p %6.4f %6.4f≤%6.4f≤%6.4f  %6.4f< <%6.4f\n"
  177.                     ,n,p,pLower,k/(double)n,pUpper,PLower,PUpper);
  178.             }
  179.         }
  180.     #endif
  181.     
  182.     double Binomial(double p,long k,long n)
  183.     /*
  184.     Returns the probability of k or more heads in n flips, where probability of
  185.     each heads is p. This identity appears in Numerical Recipes in C, page 182, and
  186.     in Abramowitz and Stegun, p. 945. Eq.26.5.24.
  187.     */
  188.     {
  189.         if(k>n)return 0.0;
  190.         if(k<=0L)return 1.0;
  191.         if(p==0.0){
  192.             if(k==0)return 1.0;
  193.             else return 0.0;
  194.         }
  195.         if(p==1.0)return 1.0;
  196.         return IncompleteBeta(p,k,n-k+1);
  197.     }
  198.     
  199.     double BinomialPdf(double p,long k,long n)
  200.     /*
  201.     Returns the probability of exactly k heads in n flips, where probability of
  202.     each heads is p. 
  203.     I'm not really sure whether this is more computationally efficient than computing
  204.     it directly, from the definition of the binomial distribution, but it is very easy
  205.     to write, and avoids the difficulty of computing the binomial coefficient
  206.     without overflow.
  207.     */
  208.     {
  209.         if(k<0L || k>n)return 0.0;
  210.         if(k==n) return pow(p,n);
  211.         return Binomial(p,k,n)-Binomial(p,k+1,n);
  212.     }
  213.     
  214.     double InverseBinomial(double P,long k,long n)
  215.     /* Returns the P-th quantile for the probability p of a heads,
  216.     given k heads in n flips. 
  217.     pUpper=InverseBinomial(0.975,k+1,n) is a 97.5% confidence upper bound on p
  218.     pLower=InverseBinomial(0.025,k,n) is a 97.5% confidence lower bound on p
  219.     Taken together these bounds form a 97.5%*97.5%=95% confidence interval: [pLower,pUpper].
  220.     */
  221.     {
  222.     #if 0
  223.         return InverseIncompleteBeta(P,k,n-k+1);
  224.     #else
  225.         /*
  226.         This simple-minded bisection routine is slow, but its answer is accurate
  227.         to within ±1e-10. 
  228.         */
  229.         double low=0.0,high=1.0,mid;
  230.         double f;
  231.         int i;
  232.         
  233.         for(i=0;i<30;i++){
  234.             mid=(low+high)/2.0;
  235.             f=Binomial(mid,k,n);
  236.             if(f>P)high=mid;
  237.             else low=mid;
  238.         }
  239.         return (low+high)/2.0;
  240.     #endif
  241.     }
  242.  
  243.     double IncompleteBeta(double x,double a,double b)
  244.     /*
  245.     The incomplete beta function Ix(a,b).
  246.     The Numerical Recipes routine assumes a>0. This work around is taken from
  247.     Abramowitze and Stegun page 944, Eq. 26.5.16.
  248.     */
  249.     {
  250.         if(a>0.0)return betai(a,b,x);    /* Numerical Recipes in C */
  251.         else return IncompleteBeta(x,a+1.0,b)
  252.             +exp(gammln(a+b)-gammln(a+1.0)-gammln(b)+a*log(x)+b*log(1.0-x));
  253.     }
  254.     
  255.     double InverseIncompleteBeta(double p,double a,double b)
  256.     /*
  257.     This simple-minded bisection routine is slow, but its answer is accurate
  258.     to within ±1e-10. 
  259.     */
  260.     {
  261.         double low=0.0,high=1.0,mid;
  262.         double f;
  263.         int i;
  264.         
  265.         for(i=0;i<30;i++){
  266.             mid=(low+high)/2.0;
  267.             f=IncompleteBeta(mid,a,b);
  268.             if(f>p)high=mid;
  269.             else low=mid;
  270.         }
  271.         return (low+high)/2.0;
  272.     }
  273.     
  274. #endif
  275.